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Abstract 



We follow a polynomial approach to analyse strong stability of linear difference 
equations with rationally independent delays. Upon application of the Hermite 
stability criterion on the discrete-time homogeneous characteristic polynomial, as- 
sessing strong stability amounts to deciding positive definiteness of a multivariate 
trigonometric polynomial matrix. This latter problem is addressed with a converg- 
ing hierarchy of linear matrix inequalities (LMIs). Numerical experiments indicate 
that certificates of strong stability can be obtained at a reasonable computational 
cost for state dimension and number of delays not exceeding 4 or 5. 
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1 Introduction 

In general, spectrum-based analysis of time-delay systems can be handled in the same way 
it is done for delay-free systems. Although the spectrum is infinite, stability is determined 
by the rightmost eigenvalues, more precisely by the sign of the spectral abscissa, the 
maximum real part of the eigenvalues. For retarded systems, the spectral abscissa is 
nonsmooth but continuous in all parameters of the system, including time delays, see [21] . 
However, it results from [121 El El E] , that, in general, it is not the case for neutral systems 
and kernel operators - the so-called associated difference equation, see also [13 [HI [19] . 
It is well-known that the spectral abscissa of the difference equation is not continuous in 
delays. Thus, arbitrarily small changes in the delay values can destroy stability. Moreover, 
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it can even happen that the number of unstable roots increases stepwise from zero to 
infinity. In order to handle this hypersensitivity of the stability of the difference equation 
with respect to delay values, the concept of strong stability was introduced by [8j. Let 
us remark that the strong stability concept has recently been generalized by [20J towar d 
difference equations with dependencies in the delays. 

As stability of its kernel operator is a necessary condition for stability of a neutral system, 
all the hypersensitivity stability issues are carried over to the stability of neutral systems. 
Thus the strong stability test should always be performed to guarantee practical stability 
of neutral systems. However, as will be shown later in the text, the strong stability test is 
rather complex. So far, a coarse numerical implementation of the test without guarantee 
or certificate has been used as a rule, see e.g. [HI [25] . Even though this brute force based 
approach works in most cases, it might fail due to approximation errors in the numerical 
scheme. As the main result of this paper we propose a more rigorous strong stability test 
that is based on a polynomial approach, relying on the numerical solution of a hierarchy 
of linear matrix inequalities (LMIs). 

In the field of time-delay systems, LMIs are usually used as stability determining criteria 
resulting from the Lyapunov time-domain approach, see e.g. [2TJ or [T6], among many 
others. 



1.1 Problem statement 



We consider a neutral system of the following form 



d ( m \ P 

— ( x{t) + H kx{t -n))= A x(t) + ^ A 3 x(t - fy) (1) 



k=l / j=l 



where x G W 1 is the state, 77. > 0, k — 1, . . .m and i?j > 0, j — 1, . . .p are the time delays. 
It is well-known, see [7i], that a necessary condition for stability of neutral system ([T]) is 
stability of the associated difference equation 



x(t) + Y,H k x(t-r k ) = 0. (2) 

k=l 

Moreover, strong stability of equation (T5]) is required, i.e. stability independent of the 
values of the delays, [21 E|- In [S] (Theorem 2.2 and Corollary 2.2), a condition for strong 
stability condition is stated as follows: 



Proposition 1 Delay difference equation |1P is strongly stable if and only if 



V, := max r a I V H k e™* < 1, (3) 
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where r a denotes the spectral radius, i.e. the maximum modulus of the eigenvalues. Fur- 
thermore, 21/70 > 1 then equation (0|) is exponentially unstable for rationally independent 
delays. 

Notice that the quantity 70 does not depend on the value of the delays, i.e. exponential 
stability locally in the delays is equivalent with exponential stability globally in the delays 

Let us remark that by homogeneity, the expression of 70 can be simplified to 

tm-l \ 
V H k e- i9k + H m ). (4) 
-1 J 

We conclude the section with some properties of the quantity 70, see [121 [IB], for more 
details. 

Properties 

1. Stability of difference equation (J2]) with rationally independent delays implies strong 
stability, and vice versa 

2. In the case of one delay (m = 1), 

70 = r CT (#i). 

3. In the case of a scalar equation (n — 1), 

m 
k=l 

4. A sufficient, but as a rule conservative, condition for strong stability is given by 

rn 
k=l 

where |.| denotes the matrix Euclidean norm, i.e. the maximum singular value. 



1.2 Computational issues 

The problem of solving (jBJ) can be formulated as an optimization task with the objective 
to find the global maximum of spectral radius over 6 e [0, 27r] m . However, in general the 
objective function r a {6) is nonconvex, i.e. it can have multiple local maxima. Besides, the 
function can be nonsmooth (e.g. at the points where the spectral radius is determined by 



-"^The m numbers t = (n, . . . , r m ) are rationally independent if and only if X^feLi n fe T fc — 0> n k € Z 
implies = 0, Vfe = 1, . . . , m. For instance, two delays t\ and Ti are rationally independent if their 
ratio is an irrational number. 
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more than either one single eigenvalue or a couple of complex conjugate eigenvalues). The 
fact that the function is nonsmooth precludes the use of standard optimization procedures. 
Instead, nonsmooth optimization methods can be used, such as gradient sampling, see 
jU [22]. However, even though these methods can handle the problem of nonsmoothness, 
they converge to local extrema as a rule. As suboptimal solutions are not sufficient (the 
global maximum of the spectral radius is needed) a brute force method has been used 
to solve the task so far, see [TS] [201 ES]- In the first step, each dimension of [0, 27r] m 
is discretized to N points. Then evaluation of ([3]) consists in solving N m times n x n 
eigenvalue problems. Hence, the overall cost of one evaluation of 70 is 0(N m n 3 ), see 
[25] . If the simplified expression 01]), the computational costs reduces to O (N m ~ l n 3 ). 
Obviously, the complexity of the computation grows considerably with the number of 
delays in the difference equation. Moreover, the risk of missing global extrema due to 
sparse or inappropriate gridding cannot be avoided. 

2 Strong stability and Hermite's condition 

Consider the characteristic polynomial 



which is homogeneous of degree n in m + 1 variables Zk, k = 0,1, . . . ,m. 

Based on (j3J), considering Zk = e? Qk , Ok £ [0, 2n], k = 1, .., m, the difference equation (j2J) is 
strongly stable if and only if the univariate polynomial 



is discrete-time stable,i.e. it has all its roots in the open unit disk. 

In order to deal with stability of this polynomial, we use a stability criterion based on the 
Hermite matrix. It is a Hermitian matrix of dimension n whose entries are quadratic in 
the coefficients of the polynomial. The Hermite matrix z\, . . . , z m — > H(z) is therefore a 
trigonometric polynomial matrix in m variables z\, . . . , z m . 

Derived by the French mathematician Charles Hermite in 1854, the Hermite matrix cri- 
terion is a symmetric version of the Routh-Hurwitz criterion for assessing stability of a 
polynomial. It says that a polynomial p(z) = p + p\z + ■ — Vp n z n has all its roots in the 
open upper half of the complex plane if and only if its Hermite matrix H(p) is positive 
definite. Note that H(p) is n-bj-n, Hermitian and quadratic in coefficients pk, so that the 
above necessary and sufficient stability condition is a quadratic matrix inequality (QMI) 
in coefficient vector p = [p p 1 ■ ■ ■ p n ] . 

The standard construction of the Hermite matrix goes through the notion of Bezoutian, 
a particular form of the resultant. A bivariate polynomial is constructed, from which a 
quadratic term is factored out, yielding a quadratic form shaped by the Hermite matrix. 
The construction is explained e.g. in [9J and references therein. See especially [15] which 



rn 




(5) 



k=l 



z ->• p(z) 
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explains that a discrete-time Hermite matrix, sometimes called Schur-Cohn or Schur- 
Fujiwara matrix, can be obtained similarly. The discrete-time Hermite matrix is also 
quadratic in the Pk, and it is positive definite if and only if polynomial p(z) has all its 
roots in the open unit disk. 

Zdenek Hurak pointed out that there is a much simpler construction of the Hermite matrix 
in the discrete-time case. The construction can be traced back to Issai Schur [23], and it 
is explained in [Tj. Entry wise formulas are also described in [3J Theorem 3.13]. Let 
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be n-by-n upper-right triangular Toeplitz matrices. Then 

H(p) = S^(p)S 1 (p) - S^(p)S 2 (p). 



Strong stability of the difference equation is hence equivalent to positive definiteness of 
the Hermite matrix of the univariate characteristic polynomial, which is a multivariate 
trigonometric polynomial matrix in Zi, . . . , z m . We express this constraint as 

H( Zl ,...,z m )yO. (6) 



3 Positivity of trigonometric polynomials 

As shown in the previous section, the key ingredient in our approach to strong stability 
of difference equation is assessing positivity of multivariate trigonometric polynomials. 
This topic has been subject to recent studies, and the recent monograph [5] is a good 
introduction focusing on signal processing applications. 

In this section we start with a scalar multivariate trigonometric polynomial, formulate 
its positivity test as a minimization problem, describe an LMI hierarchy yielding an 
asymptotically converging monotonically increasing sequence of lower bounds. We also 
describe a hierarchy of eigenvalue problems (linear algebra, much simpler computationally 
that LMI methods) to generate a hierarchy of upper bounds. 

Then we extend these results to matrix polynomials, and describe the hierarchy of LMI 
problems that must be solved to guarantee positivity of a trigonometric matrix polynomial 
at the price of solving a hierarchy of convex problems, the decision variables being entries 
of a Gram matrix yielding a sum-of-squares decomposition for the matrix polynomial. 



3.1 Minimising trigonometric polynomials 

A trigonometric polynomial has the form h(z) = h a z a where integer vector a G N n 
is a multi-index such that z a = YYl=i Z T ■> complex vector z G C™ contains indeterminates 
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such that Zi = e j9i for some 9 G [0, 27r] n , and complex numbers h a G C are coefficients. 
We use the notation z G T n to capture the constraint that each variable z k G C belongs 
to the unit disk T. 

We consider real trigonometric polynomials such that h(z) = h(z)* where the star denotes 
complex conjugation. These are such that J2 a haZ a = J2 a K* z ~ a an< ^ hence h a = h*_ a . 

Since h(z) maps T n onto R, we are interested in solving the problem 

h min = min h(z). 

zeT n 

3.2 Hierarchy of lower bounds via SDP 

In this section we construct a monotonically decreasing sequence of lower bounds on h min 
that converges asymptotically. Each bound can be computed by solving an LMI, a convex 
semidefinite programming (SDP) problem. 

First note that by definition 

h m i n = min / h(z)d/i(z) (7) 

where the minimisation is over all probability measures defined on the sigma-algebra of 
the multidisk T n , see Chapter 5 in |14j . 

Let us express polynomial h(z) as a Hermitian quadratic form 

h(z) = b* k (z)X k b k (z) (8) 

where b k (z) is a vector basis of trigonometric polynomials of degree up to k, e.g. con- 
taining monomials z a , a > 0, maXj = i ) ... |TO aij < k. Matrix is called the Gram matrix 
of polynomial h(z) in basis b k (z). Then a result of functional analysis by M. Putinar, 
transposed to trigonometric polynomials [51 Theorems 3.5 and 4.11], states that h(z) > 
if and only if there exists a finite integer d and a positive semidefinite Hermitian matrix 
X d >z such that © holds for k = d. 

As soon as k is fixed, finding a matrix X& y satisfying dHJ) can be cast into an SDP 
feasibility problem which amounts to expressing polynomial h(z) as a sum-of-squares 
(SOS) of trigonometric polynomials of degree k. 

Now defining 

h k = sup h 

s.t. h(z) — h — b* k (z)X. k b k (z) for some X. k y 

it follows that h k < h k+1 and we expect that lim^oo h k = h m m, even though a rigorous 
proof of convergence is out of the scope of this paper. 

3.3 Hierarchy of upper bounds via EVP 

In this section we show that we can construct a monotonically increasing sequence of 
upper bounds on /i min that converges asymptotically. Each bound can be computed by 
solving an eigenvalue problem (EVP) 
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In problem (J7|) let us consider that measure fi is absolutely continuous w.r.t. measure 
v, the probability measure supported uniformly on the multidisk. Let us further restrict 
the class of measures by considering that there exists a trigonometric polynomial q k (z) = 
Eo<a<i;fta za = Q*^fc( z ) °f total degree k such that (J>k(dz) = ql(dz)q k (dz)is(dz), with 
limfc^oo/ifc = // since T n is compact. Let y a = f Tn z a dv[z) denote the moment of order a 
of v. Finally, let us define 

hk = min / h(z)dfi k (z) 

as an optimisation problem over this restricted class of measures. 
With these notations 



h(z)dfi k (z) = j h(z)q* k (z)q k (z)du(z) 
h(z)qlb k (z)b* k (z)q k dp(z) 



is the same as 

where M k (hy) is called the localising matrix of order k of measure v w.r.t. polynomial 
h, see [H]. Its rows and columns are indexed by multi-indices (3 and 7 respectively, and 
its entry 7) is equal to J2 a haVa-p+y Therefore matrix M. k (hy) can be obtained from 
the moments of is, and hence it is given. It is positive definite. 

If h(z) = 1, matrix M&(y) is called the moment matrix of order k of measure v. Its entry 
7) is equal to y~p+ 7 , and hence matrix M fc (?/) is given as well. Since /i^ is a probability 
measure 

df-ik = j qtqkdvk = qlM k (y)q k = 1 

and hence _ 

h k = min qfe q* k M k (hy)q k 

s.t. q* k M k (y)q k = 1. 

It follows that h k < h k +i and lim^oo hk = h min even though I am not totally confident 
that this latter result is correct. 

Finally, given positive definite Hermitian matrices A and B, optimisation problem min,, v*Av 
s.t. v*Bv = 1 can be solved via linear algebra. Indeed, let z denote an eigenvalue of 
the pencil zB — A, and let v denote the corresponding unit eigenvector. Then vector 
v = (v*Bv)~2v is such that v*Bv = 1 and v*Av = z. Minimising this quantity then 
amounts to finding the minimum eigenvalue of pencil zB — A. 



3.4 Polynomial matrices 

The above results on scalar polynomials can be extended directly to polynomial matrices 
by considering a matrix basis instead of a vector basis to build the Hermitian matrix 
representation (jHJ). 
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In the context of our strong stability analysis problem, the core idea is then to replace the 
(typically difficult) Hermite matrix positivity condition (jH]) with a hierarchy of tractable 
SDP problems. We write 

h k = sup h 

s.t. H{z)-h = {b k {z)®L n yit k {b k {z)®I n ) (9) 

x k y 

as an LMI relaxation of order k of positivity condition fl6]). 
If Z^fc > for some k, then it implies that ([6]) is satisfied. 

If hk < f° r some k, then we cannot conclude directly, but we can try to extract from 
the dual (moment) SDP problem a certificate that indeed matrix H(z) cannot be positive 
definite, see pTj even though the trigonometric polynomial matrix case is not developed 
in this reference. If we cannot extract useful information from the dual problem, we have 
to increase the value of k and solve the next LMI in the hierarchy. 

4 Complexity 

Let M denote the size of the Gram matrix in SDP problem (Q. If we use an interior- 
point method, the worst-case complexity of one Newton iteration for an SDP problem 
in a cone of that size is 0(M 6 ). Experiments reveals that the practical complexity is 
approximately 0(M 4 ). 

The number of monomials of m variables of degree k in basis b k (z) is equal to (k + l) m . 
Polynomial p(z) has m variables and degree n so degree in (jSJ) should be such that 
2k > n. Note that we can have 2k > n since higher-degree terms may cancel in the right 
handside of equation (jSj). 

If we choose k = n/2 or k = (n + 1)/2 depending on whether n is even or not, in terms of 
complexity M = 0(ri m+1 ). The overall complexity of our SDP approach to strong stability 
analysis therefore grows exponentially in the number of delays m, and polynomially in 
the number of states n. However the exponent of this polynomial growth is quite large. 
In comparison, the gridding approach mentioned at the beginning of the paper has a 
complexity which also grows exponentially in the number of delays, but the dependence 
on the number of states is only cubic. However, contrary to the SDP approach, the 
gridding approach does not provide guarantees. 

5 Examples 

Preliminary numerical examples indicate that the EVP approach of paragraph 13.31 yields 
a sequence of bounds which converges slowly (sublinearly). This is why in this section we 
focus exclusively on the SDP approach of paragraph 13.21 

We implemented a collection of Matlab functions to manipulate trigonometric polynomi- 
als, Hermite matrices, and formulate SDP problems corresponding to positivity checks. 
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The functions are available for download^] and they provide the following functionalities: 

• sampledet .m - given a collection of matrices H^, k — 0, . . . , m, this function com- 
putes the coefficients of the multivariate polynomial p(z) = det(H + H\Z\ + ■ ■ • + 
H m z m ); it proceeds by sampling and interpolation, as described in [10] 

• trigoherm.m - computes the Hermite matrix of a homogenized multivariate poly- 
nomial; it uses the formula of ||3J Theorem 3.13] adapted to complex coefficients 

• trigohermgram.m - computes the SDP problem corresponding to the positivity test 
for a given Hermitian multivariate polynomial matrix; the SDP problem is given in 
SeDuMi's input format 

min c T x max b T y 

s.t. Ax = b s.t. z = c — A T y 
x G K z G K 

where x eR N ,y G M M , and K is the cone of positive semidefinite matrices of size 

S = y/N. 

Some instrumental functions are also provided, namely genmon . m which generates powers 
of monomials and locmon.m which locates a monomial in a Gram matrix. Besides, the 
function bfssde.m is available to evaluate ([3]) by brute force, as explained in subsection 

o 



5.1 Three states, two delays 



We adopt the illustrative example from [TS] with n = 3, m = 2, where 
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for which bfssde.m (with iV = 360) provides 70 = 0.7507 in less then 0.1 seconds under 
Matlab 7.7 on our Linux PC equipped with Intel Xeon 2.67GHz CPU with 8GB RAM. 
On Fig. [T] shows the spectral radius as a function of Q\. 

The following Matlab script assesses stability of the corresponding difference equation by 
first building the determinantal polynomial, then the corresponding Hermite matrix, then 
the SDP problem, and eventually by solving the SDP problem with SeDuMi, a primal-dual 
interior-point solver: 

H1=[0 0.2 -0.4;-0.5 0.3 0;0.2 0.7 0]; 

H2=[-0.3 -0.1 0;0 0.2 0;0.1 0.4]; 

p=sampledet ({eye (3) ,H1 ,H2}) ; °/ evaluate determinant 

p=p( : , abs (p(l , : ) ) >le-8) ; % remove small coefficients 

H=trigoherm(p) ; % compute Hermite matrix 

[A,b,c,K]=trigohermgram(H) ; % build SDP problem 

[x,y , inf o] =sedumi (A,b, c ,K) ; °/ solve SDP problem 



homepages . laas . f r/henrion/sof tware/tr igopoly . tar . gz 



9 



0.8 



CD 



X 



X 

+ 




0.4 







1 



2 



3 



4 



5 



6 







Figure 1: Spectral radius r (T (9i) for the example in Subsection 15.11 



The resulting SDP problem has size N = 2304, M = 225 and a positive semidefinite 
Gram matrix of size S = 48 is found after less than 0.1 seconds with SeDuMi 1.3. 

We can also specify the strong stability radius 70 as a second input argument to func- 
tion trigoherm. Internally, the polynomial is scaled appropriately and positivity of the 
Hermite matrix is assessed: 

H=trigoherm(p, 0.750) ; 
[A,b,c,K] =trigohermgram(H) ; 
[x,y,info]=sedumi(A,b,c,K) ; 

With the above sequence the SDP problem is found feasible. Changing the first instruction 
to 

H=trigoherm(p, 0.751) ; 

makes the resulting SDP problem infeasible, and this is certified by SeDuMi which returns 
a dual Farkas vector. As discussed at the end of paragraph 13. 41 further analysis is required 
to conclude that indeed the Hermite matrix cannot be positive definite. We leave a 
comprehensive treatment of this case for further work. 

5.2 Four states, three delays 

We consider a system with n = 4, m = 3, where 

Hl=[-0.15 0.32 0;0 -0.07 0.05; 
0.08 0.04 0;0.2 0.03 -0.13]; 
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Figure 2: Spectral radius r a (9i,8 2 ) for the example in Subsection 15.21 

H2=[-0.02 0.12 0.25;0 -0.05 0.04 0; 

0.23 -0.3; 0.19 0.28 -0.09]; 
H3=[0 -0.03 0.14; 0.01 -0.04 0; 

0.09 0.26; 0.05 -0.27 -0.06 0]; 

for which bfssde.m (with N = 360) provides 70 = 0.6028 in 4.5 seconds, see Fig. [2] with 
the distribution of the spectral radius with respect to values of 6\ and 82- 

The resulting SDP problem has size N = 250000, M = 5840 and a positive semidefmite 
Gram matrix of size 5* = 500 is found after approximately 6 minutes of CPU time, 
certifying that the spectral radius is less than one. 

5.3 Four states, four delays 

We conclude with an example with n = 4, m = 4 and the matrices 

Hl=[0.1 -0.2;pi/5 -0.1 -0.3; 

0.03 2;0 -exp(-l) 0.23] 
H2=[0 0.0456;0 -0.33 0.11 0; 
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1 0.2 0;0 -exp(-3) 0.176 0.73] 
H3=[0.1 0.65 0.42; 0.087 -pi/8 -0.1 0; 

-0.063 0.72; 0.076 0.1 -0.23] 
H4= [-0.678 -0.4; -0.0983 0; 

0.0763 0.2;-exp(-5) 0.36 0] 

for which bfssde.m (with N = 360) provides 70 = 1.7649 in more than 30 minutes. 

The resulting SDP problem has size N = 6250000, M = 52496 and a positive semidefinite 
Gram matrix of size S = 2500. This problem cannot not be solved on our computer, 
SeDuMi issues an out of memory error message. In this case, we may to try to exploit 
the problem structure (sparsity) to generate a smaller SDP problem, but this is out of 
the scope of this paper. 

6 Conclusions 

In the context of neutral time-delay systems, strong stability of difference equations is 
generally assessed numerically with a brute force gridding approach. A parallel can be 
draw with the /x-analysis approach to robustness of linear systems, see e.g. (26] where 
brute force gridding can yield misleading results and should be replaced, if possible, with 
more rigorous certificates of robustness. 

In this paper, using the Hermite stability criterion for discrete-time polynomials the prob- 
lem of assessing strong stability is reformulated as the problem of deciding positive def- 
initeness of a trigonometric matrix polynomial of size equal to the state dimension and 
number of variables equal to the number of delays. This decision problem is hard, but 
it can be approached through a converging hierarchy of tractable semidefinite program- 
ming (SDP) or linear matrix inequality (LMI) relaxations. Numerical experiments reveal 
that the approach is limited to small state dimension and a small number of delays, as 
expected. 
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